/*
 * Copyright 2010-2012 Susanta Tewari. <freecode4susant@users.sourceforge.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package bd.org.apache.commons.math.linear;

import bd.org.apache.commons.math.exception.NonSquareMatrixException;
import bd.org.apache.commons.math.util.FastMath;

import java.util.Arrays;

/**
 * Class transforming a symmetrical matrix to tridiagonal shape.
 * <p>A symmetrical m &times; m matrix A can be written as the product of three matrices:
 * A = Q &times; T &times; Q<sup>T</sup> with Q an orthogonal matrix and T a symmetrical
 * tridiagonal matrix. Both Q and T are m &times; m matrices.</p>
 * <p>This implementation only uses the upper part of the matrix, the part below the
 * diagonal is not accessed at all.</p>
 * <p>Transformation to tridiagonal shape is often not a goal by itself, but it is
 * an intermediate step in more general decomposition algorithms like {@link
 * EigenDecomposition eigen decomposition}. This class is therefore intended for internal
 * use by the library and is not public. As a consequence of this explicitly limited scope,
 * many methods directly returns references to internal arrays, not copies.</p>
 *
 * @version $Id: TriDiagonalTransformer.java 1244107 2012-02-14 16:17:55Z erans $
 * @since 2.0
 */
class TriDiagonalTransformer {

    /**
     * Householder vectors.
     */
    private final double householderVectors[][];

    /**
     * Main diagonal.
     */
    private final double[] main;

    /**
     * Secondary diagonal.
     */
    private final double[] secondary;

    /**
     * Cached value of Q.
     */
    private RealMatrix cachedQ;

    /**
     * Cached value of Qt.
     */
    private RealMatrix cachedQt;

    /**
     * Cached value of T.
     */
    private RealMatrix cachedT;

    /**
     * Build the transformation to tridiagonal shape of a symmetrical matrix.
     * <p>The specified matrix is assumed to be symmetrical without any check.
     * Only the upper triangular part of the matrix is used.</p>
     *
     * @param matrix Symmetrical matrix to transform.
     */
    public TriDiagonalTransformer(RealMatrix matrix) {

        if (!matrix.isSquare()) {

            throw new NonSquareMatrixException(matrix.getRowDimension(),
                                               matrix.getColumnDimension());
        }

        final int m = matrix.getRowDimension();

        householderVectors = matrix.getData();
        main               = new double[m];
        secondary          = new double[m - 1];
        cachedQ            = null;
        cachedQt           = null;
        cachedT            = null;

        // transform matrix
        transform();
    }

    /**
     * Returns the matrix Q of the transform.
     * <p>Q is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
     *
     * @return the Q matrix
     */
    public RealMatrix getQ() {

        if (cachedQ == null) {
            cachedQ = getQT().transpose();
        }

        return cachedQ;
    }


    /**
     * Returns the transpose of the matrix Q of the transform.
     * <p>Q is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
     *
     * @return the Q matrix
     */
    public RealMatrix getQT() {

        if (cachedQt == null) {

            final int  m   = householderVectors.length;
            double[][] qta = new double[m][m];

            // build up first part of the matrix by applying Householder transforms
            for (int k = m - 1; k >= 1; --k) {

                final double[] hK = householderVectors[k - 1];

                qta[k][k] = 1;

                if (hK[k] != 0.0) {

                    final double inv  = 1.0 / (secondary[k - 1] * hK[k]);
                    double       beta = 1.0 / secondary[k - 1];

                    qta[k][k] = 1 + beta * hK[k];

                    for (int i = k + 1; i < m; ++i) {

                        qta[k][i] = beta * hK[i];
                    }

                    for (int j = k + 1; j < m; ++j) {

                        beta = 0;

                        for (int i = k + 1; i < m; ++i) {

                            beta += qta[j][i] * hK[i];
                        }

                        beta      *= inv;
                        qta[j][k] = beta * hK[k];

                        for (int i = k + 1; i < m; ++i) {

                            qta[j][i] += beta * hK[i];
                        }
                    }
                }
            }

            qta[0][0] = 1;
            cachedQt  = MatrixUtils.createRealMatrix(qta);
        }

        // return the cached matrix
        return cachedQt;
    }


    /**
     * Returns the tridiagonal matrix T of the transform.
     *
     * @return the T matrix
     */
    public RealMatrix getT() {

        if (cachedT == null) {

            final int  m  = main.length;
            double[][] ta = new double[m][m];

            for (int i = 0; i < m; ++i) {

                ta[i][i] = main[i];

                if (i > 0) {
                    ta[i][i - 1] = secondary[i - 1];
                }

                if (i < main.length - 1) {
                    ta[i][i + 1] = secondary[i];
                }
            }

            cachedT = MatrixUtils.createRealMatrix(ta);
        }

        // return the cached matrix
        return cachedT;
    }


    /**
     * Get the Householder vectors of the transform.
     * <p>Note that since this class is only intended for internal use,
     * it returns directly a reference to its internal arrays, not a copy.</p>
     *
     * @return the main diagonal elements of the B matrix
     */
    double[][] getHouseholderVectorsRef() {

        return householderVectors;
    }


    /**
     * Get the main diagonal elements of the matrix T of the transform.
     * <p>Note that since this class is only intended for internal use,
     * it returns directly a reference to its internal arrays, not a copy.</p>
     *
     * @return the main diagonal elements of the T matrix
     */
    double[] getMainDiagonalRef() {

        return main;
    }


    /**
     * Get the secondary diagonal elements of the matrix T of the transform.
     * <p>Note that since this class is only intended for internal use,
     * it returns directly a reference to its internal arrays, not a copy.</p>
     *
     * @return the secondary diagonal elements of the T matrix
     */
    double[] getSecondaryDiagonalRef() {

        return secondary;
    }


    /**
     * Transform original matrix to tridiagonal form.
     * <p>Transformation is done using Householder transforms.</p>
     */
    private void transform() {

        final int      m = householderVectors.length;
        final double[] z = new double[m];

        for (int k = 0; k < m - 1; k++) {

            // zero-out a row and a column simultaneously
            final double[] hK = householderVectors[k];

            main[k] = hK[k];

            double xNormSqr = 0;

            for (int j = k + 1; j < m; ++j) {

                final double c = hK[j];

                xNormSqr += c * c;
            }

            final double a = (hK[k + 1] > 0)
                             ? -FastMath.sqrt(xNormSqr)
                             : FastMath.sqrt(xNormSqr);

            secondary[k] = a;

            if (a != 0.0) {

                // apply Householder transform from left and right simultaneously

                hK[k + 1] -= a;

                final double beta = -1 / (a * hK[k + 1]);

                // compute a = beta A v, where v is the Householder vector
                // this loop is written in such a way
                // 1) only the upper triangular part of the matrix is accessed
                // 2) access is cache-friendly for a matrix stored in rows
                Arrays.fill(z, k + 1, m, 0);

                for (int i = k + 1; i < m; ++i) {

                    final double[] hI  = householderVectors[i];
                    final double   hKI = hK[i];
                    double         zI  = hI[i] * hKI;

                    for (int j = i + 1; j < m; ++j) {

                        final double hIJ = hI[j];

                        zI   += hIJ * hK[j];
                        z[j] += hIJ * hKI;
                    }

                    z[i] = beta * (z[i] + zI);
                }

                // compute gamma = beta vT z / 2
                double gamma = 0;

                for (int i = k + 1; i < m; ++i) {

                    gamma += z[i] * hK[i];
                }

                gamma *= beta / 2;

                // compute z = z - gamma v
                for (int i = k + 1; i < m; ++i) {

                    z[i] -= gamma * hK[i];
                }

                // update matrix: A = A - v zT - z vT
                // only the upper triangular part of the matrix is updated
                for (int i = k + 1; i < m; ++i) {

                    final double[] hI = householderVectors[i];

                    for (int j = i; j < m; ++j) {

                        hI[j] -= hK[i] * z[j] + z[i] * hK[j];
                    }
                }
            }
        }

        main[m - 1] = householderVectors[m - 1][m - 1];
    }
}
